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In biological systems, expression dynamics that can provide fitted phenotype patterns with respect 
to a specific function have evolved through mutations. This has been observed in the evolution of 
proteins for realizing folding dynamics through which a target structure is shaped. We study this 
evolutionary process by introducing a statistical-mechanical model of interacting spins, where a 
configuration of spins and their interactions J represent a phenotype and genotype, respectively. 
The phenotype dynamics are given by a stochastic process with temperature Ts under a Hamiltonian 
with J . The evolution of J is also stochastic with temperature Tj and follows mutations introduced 
into J and selection based on a fitness defined for a configuration of a given set of target spins. 
Below a certain temperature Tg^, the interactions J that achieve the target pattern evolve, whereas 
another phase transition is observed at Tg^ < Tg^ . At low temperatures Ts < Tg^, the Hamiltonian 
exhibits a spin-glass like phase, where the dynamics toward the target pattern require long time 
steps, and the fitness often decreases drastically as a result of a single mutation to J . In the 
intermediate-temperature region, the dynamics to shape the target pattern proceed rapidly and 
are robust to mutations of J . The interactions in this region have no frustration around the target 
pattern and results in funnel-type dynamics. We propose that the ubiquity of funnel-type dynamics, 
as observed in protein folding, is a consequence of evolution subjected to thermal noise beyond a 
certain level; this also leads to mutational robustness of the fitness. 
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I. INTRODUCTION 

The function of a biological unit is generally deter- 
mined by a phenotype, which is the result of a dynamical 
process that yields a specific pattern or structure. The 
dynamics themselves are governed by a genetic sequence. 
In evolution governed by a fixed fitness condition, a phe- 
notype that gives a higher value for the fitness function 
is selected. The genes that produce such a phenotype 
are transferred to the next generation, and thus, specific 
genetic sequences are selected. However, we note, that 
the dynamical process producing a phenotype is subject 
to noise. Accordingly, the genotype-phcnotype mapping 
is generally stochastic. 

For example, genetic information determines the 
amino acid sequence for a protein, while the tertiary 
structure responsible for its functions is generated only 
by folding dynamics. By the folding process, a structure 
is formed in the protein, and this structure serves as the 
basis for the function. The genotype-phenotype map- 
ping is formed by this folding process, and this mapping 
is stochastic because of the thermal noise in the fold- 
ing process[l|. Related folding dynamics also occur in 
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t-RNA, where the influence of thermal noise on genotype- 
phenotype mapping has been intensively investigated[3- 
At a more macroscopic level, genetic information specifies 
a gene regulatory network, which determines the dynam- 
ics for the gene expression pattern, thus giving rise to 
the phenotype. This gene expression dynamics are again 
stochastic because the number of proteins in a cell is not 
necessarily very large [!]■ In fact, the stochasticity for 
gene expression of isogenic organisms has been studied 
extensively [1, 0, 0- In general, a phenotype that gives 
rise to some particular function is generated by a dynam- 
ical process that is subject to noise. Hence, phenotypes 
of isogenic individual organisms are not necessarily iden- 
tical, and therefore, they form a distribution. 

Considering that a biological function is generally a 
result of such stochastic dynamics, the dynamic process 
that shapes the function is expected to be robust under 
such stochasticity; in other words, thephenotype for the 
function will not be sensitive to noise[7|. However, such 
robustness is not a general property of dynamics. For 
example, the complex folding process of a heteropolymer 
from an arbitrary random sequence might not have such 
robustness. In this sense, the robustness could be a result 
of evolution. How can a dynamical process robust to 
noise be shaped through evolution? 

In addition to being robust to noise, a biological sys- 
tem has to remain relatively robust to mutations in the 
genetic sequence that occur through evolution; the phe- 
notype has to be rather insensitive to changes in the ge- 
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netic sequence. Are these two types of robustness cor- 
related? Does noise in the dynamic process affect the 
evolution of mutational robustness? Indeed, a possible 
relationship between robustness to noise and robustness 
to mutation has been discussed [1, [l^, [HI, fl^- H, [l3l , 
following the pioneering study by Waddingtonjla] on the 
evolution-development relationship, which is referred to 
as canalization and genetic assimilation. However, a the- 
oretical understanding of the evolution of robustness is 
still insufficient. 

Consider a dynamic process for shaping a target phe- 
notype. To have robustness to noise in such dynamics, it 
is ideal to adopt dynamics in which the target phenotypc 
is reached smoothly and globally from a variety of initial 
configurations and is maintained thereafter. In fact, the 
existence of such global attraction in the protein folding 
process was proposed as a consistency principl e by Go 
|l6l | and as "funnel" landscape by Onuchic et al. [l7l.[l8l|. 
while similar global attraction dynamics have been dis- 
covered recently in gene regulatory networks and de- 
velopmental dynamics [20|. In spite of the ubiquity of 
such funnel-like structures for phenotype dynamics, lit- 
tle is understood about how these structures are shaped 
by the evolutionary process [l|. We also address this 
question here and show that it is indeed closely related 
to the topic of robustness to noise. 

In general, it is possible to utilize a complicated model 
that agrees well with biological reality in order to answer 
the above-mentioned questions, and this will become nec- 
essary in the future. However, at the present level of un- 
derstanding, in order to understand the concepts in the 
evolution of robustness, we choose to investigate a rather 
abstract model that can be made tractable in terms of 
statistical physics, that is, a system consisting of N Ising 
spins interacting globally. Each spin can take be either 
up or down, and each configuration of spins corresponds 
to a phenotype. The fitness is given by a function of the 
configuration of some target spins, and this fitness yields 
a biological function. An equilibrium spin configuration 
is reached by a certain Hamiltonian that is determined by 
the interaction between spins. This interaction is given 
by genes and can change by mutation. By selection ac- 
cording to the fitness function, a Hamiltonian that results 
in a higher fitness is selected. Indeed, this type of model 
has been adopted by Saito et al. who utilized it in the 
study of the evolution of protein folding dynamics, where 
the spin configuration corresponds to that of residues in 
a peptide chain, and the folding dynamics are given by 
decreasing the energy in accordance with the Hamilto- 
nian. 

Even though the spin model is abstract, it can ac- 
count for the basic structures required to study the evo- 
lution of genotype and phenotype, i.e., gene develop- 
mental dynamics subject to noise — > phenotype —^ fit- 
ness. In comparison with the gene transcription network 
model utilized in the study of the evolution of robustness 
B B I the present spin model is computationally 

efficient in that Monte Carlo simulations and the meth- 



ods developed in statistical mechanics of spin systems 
can be applied to answer the above-mentioned general 
questions on evolution. In fact, we will analyze the evo- 
lution of robustness with respect to such a statistical- 
mechanical framework and define a funnel landscape in 
terms of frustration, as developed in spin-glass theory 

A shorter version of our results has already been pub- 
lished as a letter [23| . in which we propose a scenario, 
based on our numerical simulations, that the ubiquity 
of funnel-type dynamics observed in biological systems 
is a consequence of evolutional process under noise be- 
yond a certain level. In this paper, we further study the 
spin model in particular on the dependence of the re- 
sult on the number of the total spins N and the target 
spins responsible for the fitness with extensive numeri- 
cal simulations. The results suggest that there exists an 
optimal ratio of the target to the total spins to achieve 
evolution of robustness over a wide range of tempera- 
ture. In other words, some degree of redundant spins 
that do not contribute to the fitness is necessary. We 
have also carried out statistical-mechanical calculations 
of the fitness landscape, to provide an interpretation of 
the relation between the funnel dynamics and robustness 
to mutations found in our numerical simulation. These 
findings will stimulate further studies on the understand- 
ing on the evolution of robustness from statistical-physics 
viewpoints. 

This paper is organized as follows. In Sec|IIl we ex- 
plain the model setup that captures the essential features 
of the evolution. In Scc lIIIl the numerical results for the 
model are presented. First, we present the dependence of 
energy and fitness on the temperature and selection pro- 
cess. Then, we show that the evolved Hamiltonians are 
characterized by the frustration in terms of the statistical 
physics of spin systems. We classify three phases on the 
basis of the robustness of the fitness to noise and muta- 
tion, and we show that a robust system is realized at an 
intermediate temperature. The system-size dependence 
of these results is also discussed. The origin of the ro- 
bustness is studied in Sec. IIVI bv analytically estimating 
the fitness landscape by using statistical mechanics. Fi- 
nally, in Sec|Vl the conclusions and prospects for further 
development are described. 



II. MODEL SETUP 

We introduce a statistical-mechanical spin model in 
which the phenotypc and genotype are represented by 
configurations of spin variables Si and an interaction ma- 
trix Jij, respectively, with i,j = 1, • ■ • , N. The spins Si 
and Jij can take one of only two values ±1, and the inter- 
action matrix is assumed to be symmetric, i.e., Jij ~ Jji. 
A set of configurations is denoted by S for the phenotype 
and by J for the genotype. The dynamics of the pheno- 
type are given by a flip-flop update of each spin with an 
energy function, which is deflncd by the Hamiltonian for 
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a given set of genotypes, 

H{S\J) = J^j^^^j 



(1) 



We adopt the Glauber dynamics as an update rule, where 
the N spins are in contact with their own heat bath at 
temperature Ts- The Glauber dynamics, satisfying the 
detailed balance conditions, yields an equilibrium distri- 
bution for a given J: 



P{S\J,Ts 



,-PsHiS\J) 

ZsiTs) ' 



(2) 



where Zs{Ts) = Tre-'^^^'^^^-^^ and Ps = Tg\ After a re- 
laxation process, the phcnotype S follows from the equi- 
librium distribution, and it is not determined uniquely 
from the genotype J; rather, it is distributed, except 
at zero temperature. The phenotype fluctuation is com- 
puted from the Glauber dynamics, and the resulting equi- 
librium probability distribution. Thus, the degree of fluc- 
tuation is characterized by the temperature Ts- 

Next, we introduce evolutionary dynamics for the 
genotype J. The genotype is transmitted to the next 
generation with some variation, while genotypes that pro- 
duce a phenotype with higher fitness are selected. The 
time scale for genotypic change is generally much larger 
than that of the phcnotypic dynamics. We assume that 
the two time scales for the phcnotypic expression dynam- 
ics and the genotypic evolutionary dynamics are sepa- 
rated, so that the variables S are well equilibrated within 
the unit time scale of the slow variable J. Then, the fit- 
ness should be expressed by a function of the phenotype 
S that is averaged with respect to the distribution. Here, 
we define the fitness as 



(3) 



i<jet 



where (• • • ) denotes the expectation value with respect 
to the equilibrium probability distribution. The set t de- 
notes a subset of S with size t ; the members of t are 
termed as target spins. We refer to configurations such 
that all target spins are aligned in parallel as target con- 
figurations, which are assumed to give a requested appro- 
priate function. By a gauge transformation on the target 
spin and the corresponding elements of J, a choice of 
any other form of spin alignment for the fitness function, 
instead of the "ferromagnetic" configuration, yields the 
same result [2l|. The fitness can be interpreted as the 
average frequency of finding the target configurations in 
equilibrium for a given J. It should be noted that in our 
model, only the target spins contribute explicitly to the 
fitness and the remaining spins have no direct influence 
on the fitness and the selection of genes. Hence, the spin 
configuration for a given fitness has redundancy. 

The genotype dynamics are a result of mutations and 
selection, i.e., changes according to the fitness function 
following random flip-flops of genes. Hence, for a genetic 
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FIG. 1: (color online) A schematic representation of our 
model. The plane represents the genotype space, where the 
circles correspond to each configuration of J. Their fitness 
is determined through the phenotypic expression dynamics, 
which are given by Glauber-type dynamics. The landscape 
in which the height of each point corresponds to the fitness 
value is called a fitness landscape. At each generation of geno- 
type evolutional dynamics, genotypes providing higher fitness 
values are selected under the selection pressure Tj. 



dynamics, we once again adopt the Glauber dynamics by 
using the fitness instead of the Hamiltonian in the pheno- 
type dynamics, where the genotype J is in contact with 
a heat bath whose temperature Tj is different from Ts- 
In specific, the dynamics for the genotype are given by a 
stochastic Markov process with the following stationary 
distribution: 



PiJ\Ts,Tj) 



g/3j*(j|rs) 
ZjiTs,Tj)' 



(4) 



where Zj{Ts,Tj) = Tre'5-'*(-^l'^s) and /3j = Tj\ Ac- 
cording to the dynamics, genotypes are selected rather 
uniformly at high values of the temperature Tj, irre- 
spective of the fitness, whereas at low values of Tj, the 
genotypes with higher fitness values are preferentially se- 
lected. In this sense, the temperature Tj represents the 
selection pressure among mutated genotypes. 

Note that the Glauber dynamics for the genotype J is 
applied over a much longer time scale than the dynamics 
for the phenotype S; the genotype J changes only dur- 
ing the reproduction of each individual, while the spin 
dynamics proceed within a developmental time scale to 
shape the phenotype. Hence, we update after the 
spins arc updated a sufficient number of times for attain- 
ing the equilibrium configurations. In actual simulations, 
a candidate J' for the next generation is set by some flips 
of a randomly chosen from the current J, while the 
transition probability from J to J' is given by Metropo- 
lis rules, min(l,exp(/3j('I'(J) - ^'(J')))- Fig. [1] shows a 
schematic explanation of our model. 

Our model provides two landscapes: the free energy 
landscape of spins and the fltncss landscape of Js. The 
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free energy landscape is determined by a configuration of 
J and the phenotypic expression dynamics correspond to 
the relaxation process on the landscape. The fitness of J, 
'I' (J), is given by the phenotypic expression dynamics on 
the free energy landscape of spins. We call the landscape 
of fitness the fitness landscape. The evolutional dynamics 
of J correspond to a random walk to the top of the fitness 
landscape; the random walk is generated by noise whose 
intensity is given by Tj. 

A statistical-mechanical spin model similar to ours has 
been studied for protein evolution [l[ , where a genetic al- 
gorithm is used for genetic dynamics. Our model, which 
is based on two equilibrium distributions, enables us to 
conduct this study by using a Markov-chain type simula- 
tion, for which a population-based simulation developed 
by a genetic algorithm is also available. Further, analyt- 
ical tools developed in statistical mechanics are helpful 
for gaining a better understanding of the model. 



III. RESULTS 
A. Fitness and Energy 

We have carried out MC simulations of the model dis- 
cussed above and studied the dependence of the fitness 
and energy on Ts and Tj. They are given by 

^{Ts,Tj) = [^{J\Ts)]j, E{Ts,Tj) = [{H{S\J))]j, 

(5) 

respectively, where [• • • ] j denotes the average with 
respect to the equilibrium probability distribution, 
P{J ,Ts,Tj). MC sampling with temperature Ts un- 
der the Hamiltonian H and the stochastic selection pro- 
cess governed by the fitness are carried out alternately. 
In our simulations of the spin dyriamics, the exchange 
Monte Carlo simulation (EMC) [2J| is introduced to ac- 
celerate the relaxation time to equilibrium and obtain the 
equilibrium spin distribution efficiently. In this section, 
we concentrate on the analysis of the equilibrium state. 
Fig. [2] (a) and (b) show the dependence of the fitness and 
the energy on Ts and Tj, respectively, for = 15 and 
t = 3. For each generation of the genotype dynamics, the 
fitness and energy are averaged with respect to the equi- 
librium distribution over 1500 MC steps after discarding 
the first 1500 MC steps; this number of steps is sufficient 
for equilibration. The data are averaged over the last 
1000 generations. The dependence on the system and 
target size will be discussed later. For any Ts, the fitness 
decreases monotonically with Tj, but the rate of decrease 
is affected significantly by Ts- The fitness for sufficiently 
low Ts remains at a high level and decreases only slightly 
with an increase in Tj, while for a medium value of Ts, 
the fitness gradually decreases to a lower level as a func- 
tion of Tj, and eventually, for a sufficiently high value of 
Ts, it never reaches a high level. This result implies that 
the structure of the fitness landscape depends on Ts, the 
temperature at which the system has evolved. 




FIG. 2: (color online) The density plots of the fitness *( J) 
and the energy (botli are given as Eq. ((5)1) are shown in (a) 
and (b), respectively, in the Ts — Tj plane with N — 15 and 
t = 3. 



The energy function, on the other hand, shows a signif- 
icant dependence on Ts ■ While the energy is represented 
by a monotonically increasing function of Ts, for high 
Tj, it exhibits non-monotonic behavior for low Tj and 
has a minimum at Tg ~ 2.0. The configurations that 
include the target pattern at the energy minimum are 
obtained around Ts — 2.0. This non-monotonicity of 
the energy corresponds to a negative specific heat in the 
sense of standard thermodynamics. This is not possible 
in quenched spin systems with fixed J. However, the in- 
teractions J depend on the temperature Ts and Tj. It 
would be convenient to obtain an explicit formula for the 
derivative of the energy with respect to T5: 

^^^§f^^f3l{[al]j + PjCovj{{H),Coys{^,H))}, 

(6) 

where Covj{A,B) = [AB]j - [A]j[B]j, Covs{A,B) = 
{AB) - {A){B) and cr| = (H^) - (i/)2. The first term of 
Eq. ^ is the usual specific heat of the random system 
and it must be positive, and Ts-dependcncc of J comes 
from the second term, which can be negative. 

The configurations of J giving rise to the highest fit- 
ness value generally have a huge redundancy. Using a 
fluctuation induced by T5, a specific subset of the con- 
figurations of J with lower energy is selected among the 
redundant conflgurations at around T5 ~ 2.0. 



B. Frustration 

In the medium-temperature range, both a lower en- 
ergy and a higher fitness are achieved. What is the 
structure in J configurations that helps to achieve this? 
The statistical physics of spin systems tells us that a 
decrease in energy implies a decrease in the frustration 
in spin configurations. By the definition of the Hamil- 
tonian Eq. Ilj, the possible minimum energy is — C^/ 
^/N , where is the number of the spin pairs. How- 
ever, if the interaction among the three spins satisfies 
JijJjk.Jki < 0, the energy per spin cannot be minimized 
to the minimum value —C ^/y /N. Such interactions are 
said to have frustration|2ll. l2d| . Meanwhile, all the inter- 
actions satisfying JijJjkJki > do not have frustration, 
and the energy of the spin states attains the minimum 
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FIG. 3: (color online) Definition of frustration in terms of the 
parameters "l>i, $2, and $3 for the case with three target spins 
(t = 3); Si {i = 1,2,3) are target spins and the remaining 
spins are non-targets. The bold lines depict the interactions 
between target spins (s Jtt), the dashed lines depict those 
between target and non-target spins (g Jto), and the dotted 
lines depict those between non-target spins (g Joo); (a) a tri- 
angle consists of three Jtt interactions; (b) a triangle consists 
of one Jtt interaction and two Jto interactions which are con- 
nected to the same non-target spin; (c) a hexahedron that 
consists of all types of interactions. By summing up the tar- 
get spins, it is represented as a triangle with a renormalized 
interaction. 



value. However, the energetically favorable spin config- 
uration cannot be uniquely determined only by the con- 
dition JijJjkJki > 0. The spin configurations that have 
low energy should be the target configurations when both 
the decrease in energy and the increase in fitness are si- 
multaneously achieved. In our case, the target spins play 
a distinct role, and therefore, we need to quantify the 
frustration while distinguishing between target and non- 
target spins; this is in contrast to the standard spin-glass 
study. 

The interactions are divided into three categories: 
those between target spins, Jtt {{Jij \ hj & *}), those be- 
tween target and non-target spins, Jto {{Jij \ i G t, j G 
I t}), and those between non-target spins, Joo ({Jij | 
j, j 4- *})■ It can be assumed that the frustration of all 
categories decreases at intermediate Tg- To confirm this, 
we should define the conditional frustration for each cat- 
egory of spins. Fig.[3Ja) shows the minimal configuration 
consisting of the interactions in Jtt-, (b) shows that con- 
sisting of the interactions in Jtt and Jto^ and (c) shows 
that consisting of the interactions in Jto and Joo- 

We first define <I>i as the frequency of positive coupling 
among target spins, i.e.. 



*i(Ts,Tj 



t{t - 1) 



(7) 



The target configurations arc energetically preferred un- 
der ferromagnetic coupling, i.e., $1 = 1, for which no 
frustration exists among the target spins (Fig. [3]Ja)). 



Second, we define $9 as 



$2(Ts,rj 



t{t^l){N ~t) 



[EE 



J ik'Jjk 



(8) 



where t{t — 1){N — t)/2 is the total number of possible 
spins: two target spins and one non-target spin. Here, 
$2 = 1 hiiplies that no frustration exists in the interac- 
tions between target and non-target spins, and thus, the 
target configuration is at an energy minimum even when 
these interactions are included (Fig.[3Jb)). 

Lastly, as a measure of the frustration among non- 
target spins, $3 is defined as 



1 



a 



N-t 



k<lgt 



iet jet 



(9) 

where * is the total number of possible pairs of non- 
target spins. Fig. [3Jc) helps to comprehend the defi- 
nition of $3. Each non-target spin interacts with all 
the target spins, and it has t interactions that are cat- 
egorized into Jto, e.g., Ji4,J24 and J34. By summing 
up all the interactions, the frustration is computed as 
^(SLi Jii)Ji5{Yfj=i J5])- By considering all the possi- 
ble non-target spins instead of sites 4 and 5, $3 is defined 
as Eq. ([9]). If $3 is equal to 1, the frustration is not in- 
troduced by the interactions between non-target spins; 
in other words, there is no frustration globally. Hence, 
the system with $1 = $2 = $3 = 1 is in the Mattis state 
[2^ , which can be transformed to ferromagnetic interac- 
tion by a gauge transformation. 

For the interaction J, with evolution under an environ- 
ment with temperature Tg, we have computed $i,<I'2, 
and $3 by performing MC simulations. In Fig. HI we 
present contour maps of (a) ^i{Ts,Tj), (b) $2(J's, Jj), 
and (c) ^3{Ts,Tj) in the Ts - Tj plane. At sufficiently 
low Tj, the frustration parameters attain the maximum 
value 1 at the intermediate T5, where the frustrations are 
extensively eliminated, while they remain finite at low T5 
and high Ts- We define the intermediate temperature re- 
gion as Tg^ < Ts < Tg^, where the frustration parameter 
$2 equals 1. These temperatures depend on Tj, and we 
express them as T§^{Tj) and Tf{Tj). We see that for 
low Tj ( < 0.05), the phase diagram is split into three 
phases. The first one is frustrated and adapted phases 
for Ts < Tf{Tj). For Tg < Tf{Tj), all {i = 1, 2,3) 
are less than unity, and hence, the frustration remains 
for target and non-target spins. 

For Ts > Tg^{Tj), $1 equals 1, so that a target config- 
uration is embedded as an energetically favorable state 
(Fig. Hl^a)). For a finite system with finite Tj, ^j cannot 
be exactly 1. However, as long as Tj is low, the devi- 
ation of $j from 1 at the intermediate temperature is 
neghgible. In contrast to <I>i, the sum of the Jij in Jto 
and Joo fluctuates around at any Ts- Fig. [5] shows the 
averages Jto = Eietjeo and Joo = Ejgojgo 
of the summation of in Jto and Joo, respectively, at 
a low Tj (0.5 X lO^'^). As shown in Fig. [51 Jto and Joo 
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FIG. 4: (color online) Density plot of local frustrations: (a) $i, (b) $2, and (c) $3 in the Ts — Tj plane. The data are computed 
by averaging over 150 genotypes J evolved at given temperatures Tg and Tj with A*' = 15 and t = 3. 



do not deviate from at any Ts- This implies that no 
specific patterns are embedded in the spin configuration 
apart from the target spins. 

For T^^{Tj) <Ts < Tf{Tj), $2 is also equal to 1, im- 
plying that the frustration among spins is not introduced 
via interactions with a non-target spin (Fig.[4l^b)). In this 
temperature range, $3 is not always equal to 1, except 
for Ts ^ 2.0, where the Mattis state arises (Fig. IDJc)). 
When $2 = 1 and $3^1, the frustration is not com- 
pletely eliminated from the non-target spin interactions 
Joo\ this is in contrast to the Mattis state. We call such 
a J configuration "local Mattis state," and it is charac- 
terized by $1 = $2 = 1 but <t>3 7^ 1. This implies that 
the interactions J have no frustration around the target 
spins, but there is some frustration between non-target 
spin interactions. The interactions J required to form 
such a local Mattis state are obtained as a consequence 
of the evolution around Tf{Tj) < Ts < Tf{Tj) for low 
Tj, where both the fitted target configuration and lower 
energy are achieved. The T5 range in which the local 
Mattis state is stabilized becomes narrower with an in- 
crease in Tj. The phase diagram of the model is shown 
in Fig. [HI 

For Ts > Tg^{Ts), the frustration parameter <i>2 is 
less than 1, and consequently, the frustration remains, 
and the fitness ^ starts to decrease and the energy in- 
creases. Thus, neither adaptation nor energy minimiza- 
tion is achieved. The parameters $^8 should converge to 
as T5 —> 00, because for random J , the numbers of 
frustrated and non- frustrated loops are equal. In fact, at 
Ts ^ 5.0, $1 also starts to decrease. 



C. Relaxation dynamics 
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FIG. 5: (color online) Ts-dependence of the averages Jto (□) 
and IZ (O) for a fixed Tj = 0.5 x 10"^ 
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Thus far, we have computed the fitness in the equilib- 
rium state by using EMC for accelerating the relaxation 
dynamics of spins. Under standard Glauber dynamics, 
the process may require much longer time steps. Here, 
we discuss the relaxation dynamics of spins for adapted 
interactions J that are obtained from evolution under 



FIG. 6: Phase diagram of the evolved Js at TV = 15 and t = 3. 
Three types of evolved J are defined on the basis of the value 
of the fitness and $2. Their properties are summarized below 
in Table m 
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FIG. 7: (color online) (a) Relaxation dynamics of the av- 
eraged magnetization of target spins, {{mt)), averaged over 
the adapted interactions J for Ts = 10~^ (solid curve) and 
Ts — 2.0 (dashed curve). The magnetization {{mt}) is evalu- 
ated by calculating the average over 30 initial conditions for 
each J and 1000 different samples of J that are drawn from 
P{J,Ts,Tj = 10-^). 

(b)The Tg-dependence of the estimated convergent value of 
{{mt}}, m* on the right axis and the relaxation time r on the 
left. The relaxation time is estimated from the time constant 
of an exponential decay of {{mt}}. 



From an estimate of the convergent value of the target 
magnetization, m^, within a given time scale, we obtain 
the relaxation time t by fitting the estimates to the func- 
tion ((m())(s) = -|- cexp(— s/r), where s is the Monte 
Carlo step of the spin dynamics. The parameters 
and r are plotted against Ts in Fig. [Tl^b) , which shows 
that T starts to increase and decreases from 1 as Ts 
decreases below Tg^. These results imply that the inter- 
actions J whose energy landscapes are rugged, similar to 
the energy landscape of a spin-glass phase, are dominant 
for Ts < T|^, whereas those with a smooth landscape 
around the target are dominant for Tg^ <Ts< T^ . The 
latter can be interpreted as a type of funnel landscape. 
Our result supports the occurrence of transitions from 
the spin-glass phase to the funnel phase at T^^caused by 
thermal fluctuation. Note that the evolutional formation 
of funnel from rugged landscapes was also observed by 
Saito ct al. in the evolution simulations of spin systems 
for protein folding 



the condition of given Ts and Tj. The target magne- 
tization mt = \\Tlii^t^i\ is computed as a function of 
time. Note that we do not use EMC here; rather, we 
adopt standard MC to directly observe the energy land- 
scape of J adapted through evolution. We calculate the 
average of mt over J drawn from an equilibrium dis- 
tribution P{J,Ts,Tj) at Ts and Tj. Fig. [7] (a) shows 
the relaxation dynamics of {{mt)) for Ts = 10^'^(< TJ^) 
and Ts = 2.0(T|i < Ts < Tf), where ((•••)) denotes 
the average over the initial conditions randomly chosen 
and over interactions J, according to P{J,Ts,Tj). In 
the simulation, we choose a sufficient low Tj{= 10"'^) so 
that the obtained interactions have high fitness values. A 
common working temperature Tg{= 10~^) for relaxation 
is also chosen to be very low in order to examine the T5- 
dependencc of the adapted interactions J at Ts and Tj. 
By performing this simulation, the landscape properties 
of the typical J adapted at T5 are clearly determined 
from a dynamical viewpoint. As shown in Fig. [7l the 
relaxation process of {{mt)) for low T5 is much slower 
even when the working temperatures T'g are the same. 
Furthermore, the relaxation process converges to a value 
mt that is less than 1 and remains at that value for a 
long time. The deviation of {{mt)) from 1 gives the frac- 
tion of the initial conditions that fails to reach the target 
within this time span, because each mt for t = 3 is either 
1 or 1/3 depending on whether the target configuration 
is reached. Indeed, the relaxation dynamics are strongly 
dependent on the initial conditions. For some initial con- 
ditions, the spins are trapped at a local minimum, and 
as a result, the target configuration is not realized over 
a long time span. After a much longer time span, {{mt)) 
approaches 1, the equilibrium value, when the spins are 
updated under the temperature Ts, i.e., the temperature 
adopted for evolution. Such dependence on initial con- 
ditions is not observed for {{mt)) when T5 > Tg^, where 
{{mt)) approaches 1 rather quickly. 



D. Robustness to mutation 

We now examine the mutational robustness of the 
evolved genotypes in detail. The robustness to mutations 
corresponds to the stability of the fitness of J with re- 
spect to changes in the J configuration. From the geno- 
types J that are generated by P{J\Ts,Tj), mutations 
are imposed by flipping the sign of a certain fraction of 
randomly chosen matrix elements in J. The value of the 
fraction corresponds to the mutation rate /z. We evaluate 
the fitness of the mutated J'{J,fi) at Tg, i.e.. 



[^{J\J,t,)\rs)UTs,Tj) 



(10) 



where P'g = l/T'g and Z{T's,J') = Tre(-'3s-W(s| J')). The 



bracket 



S 

is almost identical to that denoted 
by [• • • ] J defined above; however, the additional subscript 
{Ts,Tj) indicates the temperatures at which the geno- 
type J evolves. If /i = and Tg = Ts, [■ ■ ■]j{Ts.Tj) is 
equal to the usual fitness defined in Eq. ([5]); [5'(J'(/i = 
0)|r^ = Ts)]j(Ts,Tj) = '^iTs,Tj). In order to distin- 
guish the mutational robustness from thermal noise, wc 
set Tj = 0.5 X 10"'^ to ensure that the fitness value with 
/I = is 1 and the working temperature is Tg = 10^^. 
The fitness averaged over 150 samples of mutated J' is 
plotted against the mutation rate ^ in Fig. 4 in (23j 
for Ts = 10-^ and Ts = 2.0. For low Ts, the fit- 
ness of mutated J exhibits a rapid decrease with an 
increase in the mutation rate, but when Ts is between 
Tg^ and Tg^, the fitness does not decrease until the mu- 
tation rate reaches a specific value. We define y.c[Ts) 
as a threshold mutation rate beyond which the fitness 
\^{J'{J, /.i)|Tg)]j(j-g drops below 1. The value has 
a plateau at Tg^ < Ts < T^ . The range of tempera- 
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Phase 


Adaptation 


Frustration 


Landscape 


Robustness 


SG 


Adapted 


Frustrated 


Rugged 


Not robust 


LMS 


Adapted 


Not frustrated 


Funnel 


Robust 


PM 


Not adapted 


Frustrated 


Rugged 


Not robust 



TABLE L Three types of J— spin-glass (SG), local Mattis 
state (LMS), and paramagnetic state (PM) — are defined by 
the value of the local frustration parameter $2 (Fig. |6]). The 
adaptation and frustration of Js have been studied in the 
previous sections, and their landscape and robustness at a 
working temperature Tg have been studied in this section. 
We have adopted a very low temperature Tg to reveal the 
energy and fitness landscape given by the interaction matrix J 
evolved at each temperature Ts- The Js evolved at Ts < Tg^ 
(these Js belong to the SG phase) have high fitness values 
and frustrations. Their adaptation is not robust to noise and 
mutation because of their rugged landscapes. The Js evolved 
at < Ts < Tf, (these Js belong to the LMS phase) 
have high fitness values and less frustrations. They are robust 
to noise and mutation, and their landscapes give funnel-type 
dynamics for spins. The fitness values of the Js evolved at 
Ts > Tc? (these Js belong to the PM phase) cannot be high, 
and the Js have frustrations. 

tures that result in mutational robustness as evolution 
proceeds agrees with the range giving rise to the local 
Mattis state, where $2 is unity. In other words, muta- 
tional robustness is realized for a set of genotypes with 
no frustration around the target spins. The evolution to 
a mutationally robust genotype J is possible only when 
the phcnotype dynamics are subjected to noise in the 
range T^^ <Ts< Tf. 

To summarize, there are three phases, local Mat- 
tis, spin glass and paramagnetic disorder ones, for the 
evolved J, as shown in Fig. [HI Each phase defined by 
the frustration parameter $2 has a distinct characteris- 
tic feature in relaxation dynamics and some robustness, 
that are summarized in Table |T1 



E. Size dependence and the existence of an optimal 
target size to obtain LMS 

To check the generality of the transition to the local 
Mattis state as well as the mutational robustness, we 
examine the model for three system sizes, N = 15, 20, 
and 30; the ratio t/N = 0.2. In Fig. [HI we compare the 
Tg-dependence of (a) the fitness, (b) the energy, and (c) 
one of the frustration parameters $2 at = 15 = 
3), TV = 20 (i = 4) and iV 30 (i 6) at a fixed 
Tj = 0.5 X lO^''. As shown in Fig. [SI the behavior of 
these quantities is qualitatively similar, and by rescaling 
the temperature Ts by a factor -x/iV, the fitness, energy, 
and $2 lines merge into a single line until the PM phase 
appears (insets of Fig. [H]) . The plateau <i>2 = 1 exists at 
all system sizes we have studied, and we show that the 
rescaled temperature Tf{N)/y/N fit together. 

The rescaling factor is y/N since the order of energy 
changes from 0{N) to 0{N^/'^) by the evolution at the 



intermediate T5 because of the existence of the local Mat- 
tis states. We have defined the Hamiltonian Eq. ^ with 
the normalization coefficient 1/\/N ; in this definition, we 
consider the average over J that would be drawn from an 
i.i.d set of Js. However, at intermediate Tg, the distri- 
bution of J deviates from the uniform distribution, and 
local Mattis states appear with high probability. As a 
result, the order of energy changes, and the temperature 
Tf (TV) is proportional to \/N . 

Next, we change the number of target spins t while 
fixing iV at 30. Fig. [Sj shows (a) the fitness, (b) the en- 
ergy, and (c) the frustration parameter $2, at = 30 
and 3 < t < 15. The threshold temperature T5 at which 
the fitness value decreases rapidly increases as t increases 
from a small value to 9 and decreases for larger t. The T5- 
dependence of the energy shows non-monotonic behavior, 
indicating the existence of the local Mattis states for all 
the values of t studied here. Similar to the threshold tem- 
perature obtained from the fitness value, the temperature 
at which the energy takes a minimum value shows non- 
monotonic behavior as t increases. Further, the range of 
the plateau with $2 = 1 is maximized at t = 7 and 8. 
Eventually, the plateau vanishes at t = 12 at least. The 
PM phase, where the adaptation is not acquired, extends 
toward the lower temperature region with an increase in 
t; the LMS phase becomes narrower as t increases. Fig.[5| 
shows that the LMS exists only up to t < 0.4A^, and for 
t > 0.4A^, a direct transition from the SG to PM phase 
occurs with an increase in T5. The LMS region in which 
robustness and high fitness can be achieved is largest at 
t ^ N/A. These findings indicate that there exists an up- 
per limit of t over which the local Mattis states cannot 
be obtained and that there exists a suitable value of t for 
stabilizing the local Mattis states for a wide range of Ts- 
The above simulation results are summarized in Fig. 1101 
where the region of local Mattis state is displayed in the 
t/N-Ts plane, by fixing N and Tj at 30 and 0.5 x 10"^^ 
respectively. 

Finally, we remark on the equilibration time of the 
evolutional dynamics of J to achieve an adapted state, 
by varying the value of t. At a fixed Tj, the time scale 
increases with the increase of t, even for a fixed size N . 
Particularly, this increase is significant in the SG phase, 
while it is moderate in LMS phase. This might remind us 
of slow relaxation of the phcnotype dynamics discussed 
in the previous section. 



IV. FITNESS LANDSCAPE 

We now explain why mutational robustness is realized 
only in the intermediate range of temperatures Tg. We 
do so by performing a statistical-mechanical calculation 
of the number of fitted states in order to obtain a sketch 
of the fitness landscape for J. We estimate the degen- 
eracy of the states with highest fitness as Ts — > 0. Let 
W{En{J)) and VF*(i?„(J)) be the number of n-th ex- 
cited states and the number of the target configurations 
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FIG. 8: (color online) The Ts-dependence of (a) the fitness, (b) the energy, and (c) the frustration parameter "I>2 at = 
15, t = 3, N = 20, t — 4 and A'' = 30, t = 6. The temperature Tj is fixed at 0.5 x 10"''. The inset of each figure shows the 
dependence of each quantity on the temperature rescaled by -\/iV. The vertical axis of the inset of (b) is also rescaled by y/N. 



(a) 1 m^-_m 




FIG. 9; (color online) The Tg-dependence of (a) the fitness, (b) the energy, and (c) $2 at A'^ = 30 and t = 3,6, 9, 12, and 15. 
The temperature Tj is fixed at 0.5 x 10"'^. The convergence value of the fitness (a) at Tg ^ cxj is 2~*'^^. 
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FIG. 10: (color online) t/N~Ts diagram of the domain where 
the local Mattis states are found at = 30 and a fixed Tj — 
0.5 X 10"'^. The points on the boundary of the domain are 
estimated by the values of Ts at which $2 drops from 1, as 
in Fig. HJc) , and the line is guide to eye. 



of n-th excited states for a given J. Then, the fitness in 



the limit as — s- is given by 



W{Eo{J)) ■ 



Z{J,Ts) 

(11) 



Accordingly, as T5 0, the highest fitness, i.e., unity, 
is achieved if and only if J satisfies the condition 
M^*(£'o(J)) = W{Eo{J)). In fact, a large number of 
J satisfy this condition besides the local Mattis state. 
For example, let us consider a Mattis state J with no 
frustration at all and introduce several changes in the 
sign of the bonds between target spins Ju, target and 
non-target spins Jto, and non-target spins Jqq. This pro- 
cedure, if applied only to the bond flips to Joo, produces 
the local Mattis states. 

First, we calculate the fitness of the Mattis states 
(<i>i = $2 = 'i'3 = 1), and we consider representative ex- 
amples of the local Mattis states ($1 = $2 = ^-.^3 7^ 1) 
and target-frustrated states that have frustration among 
the target spins and hence <I>i 7^ 1. Their fitness is 
given as the ratio of the "conditioned partition func- 
tion" Z^i{J) and the partition function Z{J), as given 
by Eq. (fTT|) . The location of the frustrations does not 
influence the partition function, but it influences Zip (J). 
To determine their fltness function, we should derive the 
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FIG. 11: (color online) (a) The Ts-dependence of the fitness 
of the Mattis state (bold line, $1 = $2 = $3 = 1), one 
of the local Mattis states (dashed line, $1 = $2 = 1,4*3 = 
0.89), and one of the target-frustrated states (dotted line, 
$1 = 0.333, $2 = 1,^3 = 0.89), plotted on the left axis. The 
dashed-dotted line represents the difference A between the 
fitness of the Mattis state and that of the target-frustrated 
state; A is plotted on the right. 

(b) The fitness of J"° adapted at Ts = 10"^ (solid line) and 
Ts = 2.0 (dashed line) with A'^ = 15 and t = 3, plotted as a 
function of the number of flipped bonds of J""" . At low 
Ts, the fitness decreases rapidly for Ua — 7 and decreases 
gradually at intermediate Ts- 

partition functions Z{J) and Z,i,{J) for three types of J. 

At first, the partition function of J with x frustrated 
interactions, Z{x,Ts) is given as 

X 

Z{x,Ts)=2c{x)Y,Cff{i,N), (12) 

where 

N—i 
71— i 

(13) 

and c{x) = e''^^~^^"'""-^2 )/\/^. This expression is vaUd 
for the case where there is at most one flipped bond at 



each site. Therefore, x should be less than N/2, and the 
expression Eq. is efficient, independent of how to 
assign the value x for Ju, J to, and Joo- 

Next, the conditioned partition function of the local 
Mattis state with x frustrated Joo interactions (denoted 
as Z^'^^^ {x,Ts)) is given as 

X 

Zk'^'ix, Ts) = 2c{x) E Cff{i, N ~ t). (14) 

i=0 

This expression is valid for the case where there is at 
most one flipped bond at each site. Accordingly, x 
should be less than [N — t)/2. The fitness of the lo- 
cal Mattis states with x frustrated interactions in Joo is 
given as Z^^^^ {x,Ts) / Z{x,Ts). Furthermore, the con- 
ditioned partition function of the target frustrated states 
that have y frustrations in Ju (implying $1 ^ 1) and x 
frustrations in Joo (denoted as Z'^^ {x,y,Ts)), is derived 
from Eq. ([13]). When a bond between a pair of target 
spins is flipped from the Mattis state and frustration is 
generated among the target spins, the energy of the tar- 
get configuration increases by 2/-\/]V; therefore, 

Z^{x,y,Ts) = Z^^'''{x,Ts) x e~W^, (15) 

where we again consider that each target spin is con- 
nected at most one flipped bond and y = I,-- - ,t/2. The 
partition functions of the target-frustrated states with 
X frustrated Joo bonds and y frustrated Jtt bonds are 
given by Z{x + y,Ts), and their fitnesses are given by 
Zl^{x,y,Ts)/Z{x + y,Ts). 

Fig. [TT] shows the Ts-dependence of the typical fit- 
ness values for Mattis, local Mattis, and target-frustrated 
states. The difference between the fitness of the target- 
frustrated state and that of the Mattis state, which is de- 
noted as A, is also plotted as a function of Ts- The value 
of fitness always approaches unity as T5 0, whereas 
such degeneracy is spht by an increase in Ts- From the 
difference in fitness between the Mattis and the target- 
frustrated states, A, the ratio of the probabilistic weight 
between these states is obtained as exp(/3jA). This sug- 
gests that fewer frustrated J states around the target, 
i.e., the local Mattis states, are preferentially selected 
only at the intermediate temperature. 

Next, we introduce frustration into J to and denote the 
constructed J as , where the superscript repre- 
sents the number of altered bonds and the subscript t 
represents the condition in which the altered bonds exist 
between the target and non-target spins, i.e., Jto- There- 
fore, the frustration parameter $2 for Jto of the state 
J"" does not equal 1 . The state J° is simply the original 
Mattis state, which has the highest fitness, whereas for 
Ua = N — 1, direct computation shows that the fitness is 
the least. Again, from a straightforward calculation, it 
can be shown that the fitness of j"° remains to be the 
highest fitness up to tIq < N/2. Hence, there is a region 
in the J-state space connected by a single point mutation 
(change of sign in a single element in J) in which the fit- 
ness retains its highest value. We refer to this region as 
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the neutral space 0, H, [2^, in the sense that a mutation 
within the region is neutral. Note that in addition to this 
construction, there is more degeneracy among the fittest 
J, as shown in Fig. [TlTa). 

We have computed how the fitness decreases as J is 
changed to leave the neutral space, for Ts 0. In 
Fig.fnTb). the fitness of J"° is plotted as a function of the 
number of altered bonds n^. By just a single point muta- 
tion, the fitness decreases suddenly to its lowest value at 
some Ha ■ This implies the existence of a clear edge in the 
neutral space. The genotype located at the edge of the 
neutral space is not robust to mutation. It is obvious that 
the Mattis state, i.e., the genotype located at the center 
of the neutral space, is robust to mutation. However, 
since the fitness of genotypes is constant throughout the 
neutral space, both the robust genotypes at the center 
of the neutral space and the non-robust genotypes at the 
edge are selected with equal weight. Then, there is no 
selection pressure to eliminate genotypes that are at the 
edge of the neutral space. Although complete degener- 
acy holds only for Ts — > 0, the above argument is valid 
for sufficiently low temperatures, and therefore, the ro- 
bustness to mutation cannot be expected to exist in the 
spin-glass phase at low Ts- This is also related to the fact 
that long time scale for equilibration in the evolutional 
dynamics is required in the spin-glass phase. As shown 
in Fig. [TlTb). the fitness value around an adapted J de- 
creases abruptly against the mutation. The number of J 
configurations on the plateau with 5* = 1 that appear at 
low Ts is roughly estimated as 2^~* because of the gauge 
transformations for N — t sites, while the total number 
of possible J configurations is 2^ . The ratio 2^~* 
is strongly suppressed as N increases. This implies that 
the evolution of J hardly finds the plateau by the local 
update. 

In contrast, at intermediate Ts, the fitness landscape 
is not neutral. For example, the fitness of J"° gradually 
decreases from its highest value with increasing tIq, as 
shown in Fig. IllT b) . There is selection pressure toward 
the genotype with Ua = 0. The genotypes with larger Ua 
have both lower robustness to mutation and lower fitness, 
but less of such genotypes are selected. Hence, the evo- 
lution toward higher fitness also induces robustness to 
mutation, as a result of the correlation between fitness 
and robustness. 

On the basis of the above argument, schematic repre- 
sentations of the fitness landscape at low and interme- 
diate Ts, together with the distribution of mutational 
robustness, are shown in Fig. ll2r a) and (b), respectively. 
The mutational robustness at an intermediate temper- 
ature observed in MC simulations, as described in the 
previous section, is thus interpreted as a consequence of 
the evolution of the fitness landscape at such tempera- 
tures. 

To confirm this schematic picture of the fitness land- 
scape, we have numerically obtained the fitness distri- 
bution of mutated Ja around the adapted J. Here we 
have computed the fitness values ^! for Js mutated from 




FIG. 12: (color online) Schematic representation of the fitness 
landscape and the probability distribution of the robustness 
to mutation estimated from the fitness landscape. 

(a) At low Ts, the genotypes that satisfy the condition 
Wii{EQ{J)) = W{Eo{J)) exist in the neutral space, where 
the fitness attains its highest value. The local Mattis states 
also exist in the neutral space, but they constitute a neg- 
ligible fraction compared to the frustrated genotypes. The 
neutral space is surrounded by non-adapted genotypes with 
sharp boundaries, as in Fig. Illf b). The genotypes around 
the edge of the neutral space that have very low robustness 
to mutation survive in an evolutionary sense, similar to the 
robust genotypes located at the center of neutral space, that 
are much fewer in number. 

(b) For Ts with T§^ <Ts < Tf, the degeneracy observed at 
low values of Ts is split, and the local Mattis states have the 
highest fitness. The fitness of the local Mattis state decreases 
continuously with increasing mutation. There is a correlation 
between fitness and robustness, i.e., the local Mattis state 
with the highest fitness is the most robust to mutation, and 
the robust genotypes are selected preferentially. 



the adapted J with the mutation rate /i = 0.1. In 
Fig. [131 we have plotted the distribution of fitness values 
at Ts = 10~3 and Ts = 2.0 for TV = 15 and t = 3. As 
shown, the fitness distribution of mutated Js at low Ts 
has two peaks &t ~ \ and 0. Hence some mutations 
are neutral, while others result in a sharp drop in the fit- 
ness to its minimal values. In contrast, the fitness of the 
mutated Js at the intermediate Ts is broadly distribute 
around 0.8. This result supports the schematic fitness 
landscape in Fig. [T^ 

At the intermediate Ts, it has been shown that the 
LMS phase vanishes at sufficiently large t as seen in 
Fig. m^c). This could be understood by the Ts and t- 
dependence of the fitness of the Mattis state. The fitness 
of the Mattis state at iV = 30 is plotted in Fig. [TH which 
is obtained from Eq. and Eq. . The fitness of the 
Mattis state as well as the LMS one at the intermediate 
temperature region decreases as t is increased. As shown 
in Fig. IHa), the fitness value decreases rapidly with Ts 
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FIG. 13: The fitness distribution of mutated J witfi mutation 
rate fj. — 0.1, wfiich are generated from the evolved genotypes 
at Ts = 10"^ and at Ts = 2.0 at iV = 15. The distribution is 
obtained by 100 types of mutated J. 



of the Mattis states - 




FIG. 14: Ts- and t-dependence of the fitness *(rs,t) of the 
Mattis state at iV = 30. 



as t increases. The drop with the increase of Ts is more 
prominent for larger t, and thus the temperature interval 
to support the LMS gets narrower with the increase of t, 
and is expected to disappear for large t. 



V. CONCLUSIONS AND DISCUSSIONS 

We have considered the evolution of a Hamiltonian sys- 
tem to generate a specific configuration for target spins 
that captures the basic features required to study the 
evolution. In this study, we adopted a Markov process, 
which is given by temperature Tj and fitness ^{J), for 
evolutional dynamics. By performing numerical simula- 
tion, we found that a specific subset of J with low energy 
and high fitness is evolved at an intermediate Ts and low 



Tj. From the statistical- mechanical viewpoint, we fo- 
cused on frustration and found that the interactions J 
that evolved at the intermediate Ts are less frustrated. 
We called these J the local Mattis states. In general, the 
less frustrated J states are robust to mutation. Hence, 
the robustness of evolving states to mutation is realized 
at intermediate temperatures Tg^ < Ts < Tg^. In other 
words, robustness to thermal noise introduces mutational 
robustness; this has also been recently discussed for gene 
regulation network models [l3|. The relevance of thermal 
noise to robust evolution is thus demonstrated. 

The mechanism by which the mutational robustness is 
achieved could be understood by a statistical-mechanical 
argument. The Ts-dependence of the fitness landscape 
was determined by explicitly calculating the fitness ^'(J) 
for a typical case. It was found that the correlation be- 
tween fitness and mutational robustness is generated at 
intermediate Ts, while it disappears at sufficiently low 
Ts- Hence, the mutationally robust interactions J are 
obtained as a result of the selection under a certain level 
of noise. The evolution of the mutational robustness at 
the intermediate Ts is confirmed by applying statistical- 
mechanical theory. 

We also found that for the interactions evolved at the 
intermediate Ts, the relaxation to equilibrium progresses 
smoothly without being stuck at metastable states. For 
protein folding, such an energy landscape was proposed 
in terms of the consistency principle b y G o [l6| and as 
a funnel-landscape by Onuchic et al. [l^. The funnel 
energy landscape has no ruggedness around the folded 
state, in contrast to the spin-glass state. On the other 
hand, far from the folded state, the ruggedness in the 
landscape remains; this is different from the global at- 
traction in the ferromagnetic or Mattis state. We note 
that the landscape in a local Mattis state that is gener- 
ated as a result of evolution at TJ^ < Ts < Tg^ is such 
a funnel landscape, as there is no frustration around the 
target spins, whereas frustration around non-target spins 
can result in ruggedness in the landscape far from the 
target configuration. Indeed, such a smooth and quick 
relaxation process is observed only for < Ts < Tg^, 
whereas the relaxation is often stuck at metastable states 
for a system evolved below Tg^ . On the basis of this cor- 
respondence between the funnel landscape and the local 
Mattis state, we expect that the funnel landscape is char- 
acterized by a state with <I>i = $2 = 1- 

Biological systems have evolved and functioned at a 
range of temperatures. Recall that for a Hamiltonian 
evolved with T5 < Tg^, a large number of time steps is 
required to reach a spin configuration having the highest 
fitness, and the number of steps seems to increase with 
the number of total spins. High fitness is not achieved for 
a low-temperature region within a biologically acceptable 
time span, i.e, within a single generation. Accordingly, 
adaptive evolution is possible only when sufficient ther- 
mal noise is present; this corresponds to T5 > Tg^. 

We have found that the funnel-like landscape evolves 
in such biologically relevant temperature regions. In fact. 
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such a landscape is commonly observed not only in pro- 
tein folding but also in gene expression dynamics 13, lOj 
and in the morphogenesis of multicellular organisms [20] ■ 
Our result implies that a landscape that allows smooth 
relaxation dynamics toward the target phenotype, such 
as the above-mentioned landscape, is realized as a con- 
sequence of dynamics that are robust to thermal noise 
as well as to mutation. We expect that this type of 
landscape that results from robustness can be found in 
general in biological evolution. It will not be restricted 
to Hamiltonian dynamics for protein folding; rather, it 
will be generally applicable in developmental dynamics. 
Indeed, recent studies on the evolution of the gene reg- 
ulation network also demonstrate that the funnel-type 
dynamics evolve at the intermediate range of noise am- 
plitude values [l3 . [ist . Our results may explain the ubiq- 
uity of such funnel-like dynamics in evolved biological 
systems. 

Although the frustration measure in the present paper 
is not directly applied to the gene regulation network, 
the transition to the robust developmental landscape is 
common. It will be interesting to study the similarities 
and differences in the transitions in the spin Hamiltonian 
model presented here and the dissipative gene expression 
dynamical model. 

It is also interesting that there is an optimal number 
of target spins for achieving the local Mattis state with 
a funnel energy landscape over a wide range of temper- 
atures. The number of non-target spins N ~ t controls 
the redundancy of the system. If t is too small, the tar- 
get configurations will be more easily perturbed by non- 
target spins, and fitted states will be less robust to the 
thermal noise. On the other hand, if it is too large, such 
J configurations with high fitness are too much limited 
in the whole configurations, and hence it will not be eas- 
ily accessed or may be destabilized by mutation. The 
existence of an appropriate number of redundant spins 
is important to achieve robustness. Indeed, in proteins, 
amino-acid residues that arc responsible for a function 
are limited and their number is typically smaller than 
the number of other proteins. Possible links between re- 
dundancy and the evolvability are pointed out in , but 
they still wait to be established quantitatively by theoret- 
ical analysis. The present demonstration of the optimal 
fraction of target elements for realizing robust functions 
will be important in this context. 

In this study, we observed transitions at and TJ^. 
The phase below Tg^ corresponds to the spin-glass phase 
and that above Tg^ corresponds to the paramagnetic 
phase in the context of statistical physics; in contrast, 
the local Mattis phase between the two phases, which 
could correspond to the funnel landscape, is a novel dis- 
covery in this study. Since a framework of statistical 
physics of spin systems has been adopted in the present 
study, the theoretical concepts developed therein, such 
as replica symmetry breaking, may be applicable in the 
context developed here to understand this transition. In 
particular, our model in this study is a variant of spin 



Hamiltonian systems with two temperatures: one for spin 
and the other for interactions J . Theoretical analysis of 
such systems [23, will be relevant for the analysis of 
the local Mattis state and mutational robustness that we 
have discussed in this paper. 
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